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Abstract. The non-equilibrium phase transition in models for epidemic spreading 
with long-range infections in combination with incubation times is investigated by 
field-theoretical and numerical methods. In this class of models the infection is 
assumed to spread isotropically over long distances r whose probability distribution 
decays algebraically as P(r) ~ r~ d ~ cr , where d is the spatial dimension. Moreover, a 
freshly infected individual can infect other individuals only after a certain incubation 
time At, which is distributed probabilistically as P(At) ~ (At) _1_K . Tuning the 
balance between spreading and spontaneous recovery one observes a continuous phase 
transition from a fluctuating active phase into an absorbing phase, where the infection 
becomes extinct. Depending on the parameters a and k this transition between 
spreading and extinction is characterized by continuously varying critical exponents, 
extending from a mean-field regime to a phase described by the universality class of 
directed percolation. Specifying the phase diagram in terms of a and k we compute 
the critical exponents in the vicinity of the upper critical dimension d c = cr(3 — by 
a field-theoretic renormalization group calculation and verify the results in one spatial 
dimension by extensive numerical simulations. 
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1. Introduction 

Epidemic spreading within a population of spatially distributed individuals is a 
complex phenomenon which depends on many parameters such as the density and 
the susceptibility of individuals, the specific transport mechanism, incubation times 
and memory effects like immunisation. Depending on the interplay of these factors 
an infectious disease may either spread over the population or disappear after some 
time. In order to understand the transition between survival and extinction of the 
pathogen various models with different degrees of simplification have been discussed in 
the literature [1-5] . The interest of physicists and mathematicians in models of epidemic 
spreading is mainly motivated by the observation that in many cases the transition is 
continuous and characterized by universal properties determined by the symmetries of 
the model which are independent of microscopic details. This allows one to classify phase 
transitions in epidemic models into universality classes characterized by certain critical 
exponents and universal scaling functions [6-8]. Typically these classes correspond to 
different underlying field theories that describe the long-range properties of the dynamics 
at the transition. 

Directed percolation (DP), also known as contact process, is probably the simplest 
and best-studied universality class of transitions in epidemic models [9, 10]. This 
class is characterized by three independent critical exponents /3,u±, and vn. In the 
stationary active phase these exponents are related to the density of infected individuals 
P ~ (p~ Pc) i the spatial correlation length £j_ ~ (p— p c ) ~ u± , and the temporal correlation 
length £|| ~ (jp — p c )~ u ^, where p is the infection probability and p c its critical threshold. 
Directed percolation corresponds to a specific type of field theory which made its first 
appearance in the context of high-energy physics [11-13]. This observation led Janssen 
and Grassberger to their celebrated conjecture that any model which complies with the 
symmetry properties of this field theory should belong to the DP class [14,15]. One of 
these conditions is that the infections have to be local in space and time. 

In order to study more realistic models which are still simple enough to be 
treated analytically but sufficiently complex to exhibit non-DP behavior, two important 
generalisations have been investigated. The first one is an epidemic process with long- 
range spreading, where the infectious disease is transported in random directions over 
long distances [2]. Motivated by empirical data the spreading distances r in these 
generalized models are assumed to be distributed as a power law 

P(r)~r- d " CT , (<7>0) (1) 

where a is a control exponent. In the literature such algebraically distributed long-range 
displacements are known as Levy flights [16] and have been studied extensively e.g. in 
the context of anomalous diffusion [24]. In the present context of epidemic spreading 
long-range flights do not destroy the transition, instead they change the critical behavior 
provided that a is sufficiently small. More specifically, it was observed both numerically 
and in mean field approximations that the critical exponents change continuously 
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with a [17-19]. As a major breakthrough, Janssen et al introduced a renormalizable 
field theory for epidemic spreading transitions with spatial Levy flights [20], computing 
the critical exponents to one-loop order. Because of an additional scaling relation only 
two of the three exponents were found to be independent. These results were confirmed 
numerically by Monte Carlo simulations [21]. As a further step towards more realistic 
models, Brockmann and Geisel [22] recently studied epidemic spreading by Levy flights 
on a inhomogeneous support in the supercritical phase. 

The second generalization introduces a similar long-range mechanism in temporal 
direction. Such 'temporal' Levy flights can be interpreted as incubation times At 
before an infected individual can subsequently infect others. As in the first case, these 
incubation times are assumed to be algebraically distributed as 

P(At) ~ At- X - K , (k > 0) (2) 

where k is a control exponent. However, unlike spatial Levy flights, which take place 
equally distributed in all directions, such temporal Levy flights have to be directed 
forward in time. A field-theoretic renormalization group calculation by Jimenez- 
Dalmaroni and Cardy [23] allows the critical exponents to be computed to one-loop 
order. Again only two of the three exponents turn out to be independent because of an 
additional scaling relation. 

In the present paper we investigate the mixed case of epidemic spreading by spatial 
Levy flights combined with algebraically distributed incubation times. As we will show, 
the corresponding field theory renders two additional scaling relations, hence only one 
of the three exponents is independent. We present a phase diagram in terms of the 
control exponents o and k, specifying the mean-field phase, a DP phase, two phases 
corresponding to the previously studied cases, and a novel fluctuation-dominated phase 
describing the mixed case in which the critical exponents are computed by a field- 
theoretic renormalization group calculation to one-loop order. Investigating several cross 
sections of the phase diagram by high-precision Monte Carlo simulations we demonstrate 
that the field-theoretic predictions are in fair agreement with numerical results. 

2. Modelling long-range spreading and incubation times 

Ordinary DP with short-range interactions is described by the Langevin equation [14] 

rd t p = ap- bp 2 + DV 2 p + £(x, t) , (3) 

where p(x, t) is a coarse-grained density of infected individuals and r, a, b, D are certain 
coefficients (r is introduced for later convenience and may be set to 1). The Gaussian 
noise £(x, t) describes density fluctuations. Its intensity varies with the density of active 
sites, as described by the correlation function 

(£(x, t)£(x', f )) = 2cp(x, t) 5 d (* - x!)6(t - . (4) 
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Following Refs. [20,21,23], the generalized case of algebraically distributed long-range 
infections and incubation times is described by the Langevin equation 

(rd t + fd?)p = ap- bp 2 + {DV 2 + DV a )p + £(x, t) (5) 

with the same noise as in Eq. (JHJ). Here d£ and V 7 are linear non-local operators acting 
on p(x, t) by 

^p(x,t) = -^y o dt'r 1 - K [p(x,t)-p(x,t-t')] (0<*<1) (6) 

and 

V CT p(x, t) = jf~ ± j d d x> W\- d ~ a [p(x + x, t) - p(x, t)} (0 < a < 2) (7) 

where A/j| = — T(— k) and A/j_ = — 2 cos(7T<r/2)r(— a) are appropriate normalization 
constants. Acting on a plane wave these operators fulfill the properties 

dfe™ = {iu) K e luJt , (8) 
V CT e ik ' x = - |k| CT e ik ' x (9) 

so that they can be interpreted as fractional derivatives [24] . As an important difference 
the temporal operator is directed in time and thus brings down the nonsymmetric factor 
(— iu) K in front of the exponential, while the undirected spatial operator renders the 
factor — |k| CT which is symmetric under reflections. We note that the operator V 2 is still 
non-local and hence differs from the ordinary Laplacian V 2 although both operators act 
in the same way on a plane wave. Similarly, d] is non-local and differs from the ordinary 
local derivative d t . 

Since the integral kernels in Eqs. (jZJ) and (JHJ) would diverge for t' — > and x' — > 
one has to introduce lower cutoffs for the spreading distances in space and time. These 
cutoffs induce short-range components which are accounted for by retaining the local 
operators dt and V 2 in Eq. (JHJ). Even if these terms were not initially included, they 
would be generated under renormalization of the theory. 

A dimensional analysis of Eq. (J3J) immediately yields the mean field exponents 

oMF -i „MF -1 MF -1 / ln \ 

P = A , Z/j_ = (J , Vn = K , {10) 

and the upper critical dimension 

d c = 3a-- (11) 

below which fluctuation effects are expected to become relevant. These exponents vary 
continuously with a and k. As expected, for a — > 2 and k — > 1 they cross over smoothly 
to the mean field exponents of DP in d > 4 spatial dimensions. 
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Figure 1. Feynman diagrams of directed percolation. The figure shows the 
one-loop expansion of a) the propagator and b)-c) the cubic vertices. The thin 
lines represent the free propagator Go(k, u) directed in time according to the 
arrow while the small circles denote bare cubic vertices ±g. 

3. Field-theoretic renormalization group calculation 

3.1. Field-theoretic action 

Starting point of a field-theoretic renormalization group calculation is an effective 
action S which is defined as the partition sum over all realizations of the density field 
ip('x,t) := p(x, t) and the noise field £(x, t), constrained to solutions of the Langevin 
equation (jSJ). Integrating out the noise by introducing a Martin-Siggia- Rosen response 
field ipfeit) an d symmetrizing the coefficients of the cubic terms by rescaling the field 
ip and z/5 one arrives at the field theoretic action 

Sty,$]=S W,$]+S ist W,$] (12) 

consisting of a free contribution 

SofoM] = J d d xdti>(rd t + fd? -a- DV 2 - DV a )ip (13) 

and an interaction part 

S in thM = 9 f d d xdt (^ 2 -V5» (14) 



with the coupling constant g = ybc. 

It is important note that the additional long-range operators <9f and V CT appear 
only in the free part of the action. Therefore, the structure of the Feynman diagrams 
in a loop expansion is exactly the same as in DP (see Fig. QJ, the only difference being 
that the free propagator in momentum space G® p (k, to) = (Dk 2 — a — iru)^ 1 is replaced 
by the generalized counterpart 



Dk 2 + Dk a -a-iruj - f(iu) 



(15) 



Apart from this modification the loop integrals therefore have exactly the same structure 
as in DP. In particular, the so-called rapidity reversal symmetry 

V>(k,w) -»• -^(k,-w) i/>(k,u;) -»• -ip(k,-u) (16) 
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is still valid. As a consequence the fields ip an d if> scale identically and the cubic vertices 
shown in Fig. ^3-^: renormalize exactly in the same way. 

Wilson's momentum shell renormalization procedure consists of two steps. At first 
the action is rescaled by 

x -> 6x, t -> b z t, ip -> 6 X V, $ -> & X V> (17) 

where b = 1 — i with < £ <C 1 is a dilatation factor while z = v\\/v±_ and 
X = — /3 / u± are critical exponents. Introducing an upper cutoff A in momentum space 
this scaling transformation inflates the cutoff by A — > A/6. In the second step of Wilson's 
renormalization procedure, the newly arising short wavelengths within the momentum 
shell A < |k| < A/6 are integrated out in the loop expansion and expanded to lowest 
order in k and u, renormalizing the coefficients of the propagator and the vertices. In 
the present case, shell integration to one-loop order renormalizes the propagator and 
the coupling constant by 

G-\k,u) ->CT 1 (k,a;)+€.L(k,a;), g^g-tV (18) 

where 

L(k -» = 292 L +k '^ + ^ - - ^ (l9) 

V =8 9 3 y kA |^Gg(k',a,')G„(-k',-u,'). (20) 

are the shell integrals corresponding to the one-loop diagrams shown in Fig^ At the 
fixed point of rescaling combined with shell integration the system is assumed to be 
critical and the two exponents x an d z can be calculated. 



3.2. Directed percolation 



Before calculating the critical exponents for generalized epidemic spreading with spatio- 
temporal long-range interactions let us briefly recall Wilson's renormalization group 
approach applied to DP [25]. Setting f = D = the action involves four coefficients 
g, a, r, and D. Under combined action of rescaling and shell integration these coefficients 
change as 



dilng = 3x + d + z — L Q 
dl In a = 2x + d + z — L a 
dilnr =2 X + d-L T 
d£\nD = 2 X + d + z-2- 
with the loop corrections 



Lg , 



V 

9 



_ L(0,0) _ i dL(0,u) 

— 7 L>t — T. 

a t ouj 



L D = 



LJ = 



1 d 2 L(k,0) 
2D 



dk z 



(21) 
(22) 
(23) 
(24) 

(25) 



k=0 
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In the case of DP the shell integrals can be computed, leading to 

2g*A d K d g*A d K d L g L g 

Lg ~ r(ZM 2 -a) 2 ' a ~ ar{Dh?-ay ~ T' ~~ ~8~' W 

where K d = ^T^^fuHj is the surface of a d- dimensional unit sphere divided by (2n) d . 
Inserting these results the flow equations have the fixed point 

L* = 2 + — , L* = — , x = -2 + — , z = 2-— , (27) 
12 9 3 12 12 

where e = <i c — d denotes the distance from the upper critical dimension d c = 4. Assuming 
D and r to be stationary and expanding these relations to first order in e this fixed point 
corresponds to 

3K d 6 

In order to compute the third critical exponent one has to consider the flow equations 
for a and g in the vicinity of the fixed point. To this end one has to linearize the flow 
equations and to compute the eigenvalues of the matrix 

n, I d g9(3x + d + z - L g ) d g a(2x + d + z - L a ) \ 



\ ^(3% + d + z - L g ) d a a(2x + d + z - L a 



g=g*,a=a* 



evaluated at the fixed point and to select the positve eigenvalue 2 — e/4, which is equal 
to uj_ . To first order in e we therefore arrive at the well-known result for DP 

/3=l~ + 0(e 2 ), ^|| = l + ^ + 0(e 2 ), ^ = \ + ^ + 0(e 2 ). (30) 



3.3. Spatio-temporal long range interactions 

After summarizing the known case of DP let us now discuss the generalized model with 
spatio-temporal long-range interactions. As discussed above, by introducing long-range 
flights and incubation times we modify solely the free propagator while the structure 
of the loop integrals (fTT)|) and f!2()l) remains unchanged. However, expanding these 
integrals as a Taylor series around k = ui = one always obtains integral powers of 
k and uj while fractional terms of the form k a or (iu) K are not generated. This means 
that the fractional operators d t and V CT do not renormalize themselvses, instead they 
renormalize their short-range counterparts d t and V 2 . Assuming this observation to 
hold to any order of perturbation theory, the RG flow equations for the coefficients of 
the long-range operators are determined exclusively by the rescaling-part 

<%lnf = 2x + d + z — kz , (31) 
d e \nD = 2 X + d + z-a. (32) 

At the fixed point the right-hand sides of these equations are zero, implying two exact 
scaling relations, namely, the temporal relation 2x = o~ — d — z and the spatial relation 
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2x = i K ~ I) 2 — d. In terms of the standard exponents (3, v±, uv, these scaling relations 
read 

2p+(a-<t)vj_-v\\ =0, (33) 

2/3-dv± + (K-l)u\\ = 0. (34) 

Hence two of the three critical exponents are already determined by trivial 
renormalization of the long-range operators. The first one is the exponent for the 
temporal decay of the density 

. (3 dn + a-Ka . . 

o = — = . (35) 

u\\ 2a K ' 

The second one is the dynamical exponent z, which relates spatial and temporal 
correlation lengths: 

z=-. (36) 

K 

We note that the same result for z was obtained in the case of pure anomalous diffusion in 
space and time [24] . The combined effect of spatial and temporal long-range interaction 
apparently predominates the critical dynamics to such an extent that the dynamical 
exponent z locks onto the strict ratio ct/k. 

To obtain the third exponent one has to consider the remaining flow equations 

de In g = 3% + d + z — L g (37) 

d e \na = 2 X + d + z- L a (38) 

in the vicinity of the fixed point 

3 K *-d K -a = e 
9 2k 2' v ; 

where e = d c — d denotes the distance from the upper critical dimension 

d e = 3cr - - . (40) 
K 

To compute the matrix M in Eq. (|29|) we note that 

--3L g . (41) 



d(aL a ) _ L g d(aL a ) _ 2aL a d(gL g ) 



da 2 dg g dg 

Here the first equation relies on the identity 

^[G (kV)G (-k\-cy)] = G (kV)G 2 (-k',-u/) + (42) 

relating the integrands of Eqs. (|19|) and ()20j). Again the off-diagonal entries of the 
matrix vanish to first order in e, hence to one-loop order the positive eigenvalue of M 
is z/J 1 = —{2x + d + z + d a (aL a )) = a — L*/2 = a — e/4. Summarizing these results, we 
therefore arrive at the critical exponents to one-loop order 

0=l-£ + Otf), "ll4 + 4^ + °< f2 >- "x = ! + £+<*»).(«) 
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Figure 2. Phase diagram for epidemic spreading with spatio-temporal 
long-range interactions in one spatial dimension. The five phases and their 
boundaries are explained in the text. The circles along the grey lines indicate 
the points where the numerical simulations were carried out. 



3.4- Phase diagram 

The field-theoretic results suggest a generic phase diagram in terms of the parameters 
a and n. Fig. |21 shows the phase diagram in the case of one spatial dimension d — 1, 
for which the numerical simulations are performed (see Sect. |3J). It comprises seven 
different phases, namely: 

(i) a mean field phase (MF) for small a and k which is governed both by Levy flights 
and incubation times. In this phase fluctuation effects are irrelevant and the expo- 
nents are given by (/3, v\\, u±) = (1, re -1 , cr _1 ), see Eq. (jTHjl . 

(ii) a mean field phase for k > 1 governed solely by Levy flights (MFL). In this regime 
waiting times described by the operator d£ are irrelevant compared to dt and the 
exponents are given by (ft, v\\, u±) = (1, 1, o" -1 ). 

(iii) a mean field phase for o > 2 dominated by incubation times (MFI). Here spatial 
Levy nights associated with the operator V CT become irrelevant and the exponents 
are given by (/?, v\\,v±) = (1, k^ 1 , 1/2). 
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(iv) a directed percolation regime (DP) in the limit of short-range interactions (large a 
and k) with the exponents /3 DP ~ 0.2765, — 1.097, and z/° p ~ 1.734. 

(v) a non-trivial fluctuating phase dominated by spatial Levy flights (L), where incu- 
bation times are still irrelevant so that the results of Refs. [20,21] apply. 

(vi) a non-trivial fluctuating phase dominated by incubation times (I), where spatial 
Levy flights are irrelevant so that the results of Ref. [23] apply. 

(vii) a non-trivial fluctuating mixed phase in the center (LI), on which we focus in the 
present work. In this phase spatial Levy flights and incubation times are both 
relevant and the field-theoretic results derived in the previous subsection apply (see 
Eq. (USD). 

Assuming that the critical exponents between the phases change continuously we can 
compute the phase boundaries by means of the scaling relations (jHH|) and For 
example, the spatial scaling relation ()34j) holds inside the phases L and LI while it is 
violated outside. Hence, plugging in the numerically known exponents of the adjacent 
phases we can compute phase boundaries in terms of a and k. For example, inserting 
DP exponents into Eq. (|3*l~j) renders a vertical line at o = 1 + z DP — 2/5 DP /z/^ p ~ 2.077 
as a phase boundary. Similarly, inserting the mean field exponents (jlO)) into Eqs. (|3*3|) 
and ()34|) gives the exact relation 

«=o-^T' (44) 
3a — l 

which defines the curved phase boundary between the mixed and the mean field 
phase (line A-B in Fig. |2J). Finally, the remainig two curved phase boundaries of the 
mixed phase (lines A-C and B-C) are obtained by inserting the numerically estimated 
exponents for mere Levy flights [21] or mere incubation times [23] into the corresponding 
scaling relation. 

In dimensions d = 2, 3 the area of the mixed phase (LI) is expected to shrink but 
the phase diagram still has the same qualitative structure. In d > 4 dimensions we 
expect the phases L,W, and LW to disappear so that the phase diagram shows only four 
different mean field phases separated by crossed straight lines at cr = 2 and k — 1. 

3. 5. Critical behavior in the vicinity of d c = 1 

The predictive power of field-theoretic loop expansions is restricted to a small regime in 
the vicinity of the upper critical dimension. On the other hand numerical simulations of 
absorbing phase transitions are often carried out in one spatial dimension, i.e., usually 
far away from the regime where the field theory is valid. In the present model, however, 
the upper critical dimension depends continuously on the parameters a and k and thus 
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we can tune them in such a way that d c is close to 1. Therefore, we can quantitatively 
verify the field-theoretic predictions in 1+1-dimensional numerical simulation. 

As an example we consider first-order corrections of the critical exponents when 
moving away from the mean field phase MF into the mixed phase LI along a horizontal 
line. For given k the line where d c = 1 is given by 

*« = s£t- (45) 

Let us now choose a — a c + ji slightly above a c , where < /i C 1. To first order in /j, 
the critical exponents are then given by 



(1 - 3k) 



2 



/? = !- , , H + OQi 2 ), (46) 



4k 2 

„3-i-^ + W . (48, 
Similarly for given a we have d c — 1 if k is equal to 

k c = °" . (49) 

Choosing again k = k c + r] slightly above k c , where < rj <C 1 is the vertical distance 
from the phase boundary, the critical exponents are given to first order in 77 by 

P = 1 " + 0( V 2 ) , (50) 

n 1 (a + 1)(1 - 3a) 2 2N 
^11 = 3 - - - ^ ^_ — + 0(^ 2 ) , (51) 

^i + il^. + CV). (52) 



4. Numerical analysis 

4-1. Seed simulations 

Modelling long-range interactions on a finite lattice almost inevitably evokes significant 
finite-size effects. Simulations studying the decay of a fully occupied state usually 
operate on a confined lattice, often with periodic boundaries. However, this approach 
requires huge lattice-sizes in order to minimize finite-size effects caused by long-range 
interactions, especially if a is small. As observed in Ref. [21] this is already an issue 
when interactions are still local in time. Therefore we prefer so-called single-seed 
simulations [26] , where a cluster is generated by a single active site placed at the origin. 
Instead of a complete lattice, we now only store the individual coordinates of active 
sites. Thus, the lattice-size of the model is merely given by the machine limits for 
integer numbers (±2 63 ), virtually eliminating finite-size effects. 

In seed simulations one usually measures three time-dependent quantities, namely 
the survival probability P s (t), the current number of active sites N(t), and their current 



Epidemic spreading with long-range infections and incubation times 



12 



mean square spreading from the origin R 2 (t). In most models with absorbing states, the 
survival probability scales as P s (t) ~ t~ s , where S — (3/v\\. In the present model with 
incubation times, however, the survival probability is no longer well defined since there 
may be intermediate time intervals where the system is not active (see Ref. [23]). We 
therefore suggest to measure the autocorrelation function P a (t) which is defined as the 
probability of the site at x = being active at time t and which remains well defined 
despite of incubation times. The autocorrelation function represents the full propagator 
of the field theory and is expected to scale as 

Pa(t) ~ t- 25 (53) 

with 5 = (diK + cr — na)/2a inside the mixed phase. 

The number of presently active sites N(t) is measured as usual. This quantity is 
expected to scale as 

N(t) ~ t 9 , (54) 

where 9 = (dv± — 2j3)/v\\. Inside regions I and LI, where incubation times are relevant, 
the additional scaling relation (jHlj) implies the simple exact expression 

9 = k - 1 . (55) 

Finally, the mean square spreading from the origin is expected to scale as 

R 2 (t) ~ t 2/z , (56) 

with z = o~/k inside the LI regime. However, the arithmetic average of R 2 {t) is 
problematic as it immediately diverges for o < 2. As suggested in [21] we therefore 
consider the geometric average instead, which is finite and turns out to comply with Ij5fi)l . 

Note that because of 9 = —25+d/ z one of the three exponents is redundant. In fact, 
it would suffice to measure only two quantities (accounting for two scaling relations), 
but we deem it useful to have an additional numerical validation. 

4-2. Lattice model and its algorithmic implementation 

As a generalization of directed bond percolation, our simulation features a 1+1- 
dimensional lattice with integer positions x G Z (see FigEJ). Starting with a single 
active seed located at x — 0, the infection spreads as follows: For each active site at 
a given time, two random flight distances are generated, one directed in positive, the 
other in negative x-direction. The flight distances are distributed as 

n/ x f or' 1 ' if r > 1 

P(r) = < , (57) 

v ' 1 otherwise v ' 

This distribution can easily be generated by drawing a random number q from a uniform 
distribution between and 1 and setting the flight distance r = q^ 1 ^. Note that this 
implies a lower cutoff equal to the minimal flight distance of one lattice spacing. 
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Figure 3. Lattice model used for the simulation. An active seed at the 
origin (red) infects two target sites (green) with probability p, spreading over a 
randomly generated spatial distance r after an incubation time At (see text). 



Likewise, two random incubation times are generated, using a similar probability 
distribution determined by k: 

P(A t ) = ( KAr "" if f (58) 
v 1 | otherwise v 1 

In both cases, fractional digits are truncated in order to retrieve integer coordinates. 
Finally the two target sites are independently infected with the probability < p < 1. 

The spatial and temporal coordinates of active sites are stored in a dynamically 
generated list. In fact, we use a heap where we have immediate access to the next 
active site on the time scale to be processed. Active sites with the same time coordinate 
additionally use their spatial coordinate to maintain strict order. This allows the active 
sites to be processed in proper sequence and hence tracking down multiple activation 
easily (by definition each site can only be active once at any given time). 



4-3. Simulation procedure 

Essentially, the model features three control parameters, namely a, k and p. To begin 
with, we pick out a point in our phase diagram, specified by certain a and k. Then 
we determine the critical threshold p c for this pair of control exponents. Averaging 
over numerous independent realizations of the process, we measure the quantities P a {t), 
N(t) and R 2 {t) as introduced above. In order to determine the threshold probability 
p c , we first adjust p such that we obtain asymptotically straight lines in a double 
logarithmic plot. As usual, deviations from criticality lead to a curvature of all the 
observed quantities in a double logarithmic plot, N(t) seeming to be the most sensitive 
to this. We observe that the scaling relation 9 — k — 1 holds accurately inside the 
Li-phase (Fig. |2J). Therefore we then use this scaling relation as a fit for N(t) in order 
to improve the estimates of p c . 
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Figure 4. Estimates of the critical exponents along the horizontal sections 
shown in Fig.[2]for k = 1 (upper panel) and k = 0.5 (lower panel). Theoretical 
predictions are shown as dashed lines (see text). 



4-4- Numerical results 

In order to obtain an overall impression, we opt for four different cross sections through 
the phase diagram, indicated as grey lines in Fig. |21 Having determined the critical 
threshold for appropriate combinations of a and k along these cross sections, we measure 
P a (t), N(t) and R 2 (t) at criticality. In a double logarithmic plot the asymptotic slopes 
of these quantities yield the exponents —25, 9 and 2/z respectively. 
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In order to verify the field-theoretic one-loop expansion, it is also necessary to 
study the off-critical regime. Choosing p slightly below the critical threshold we again 
measure the average number of particles N(t). According to standard scaling arguments 
this quantity should scale as 



where / is a universal scaling function. Therefore, plotting N(t)t 9 versus t(p c — p) y w 
for different values of p we can determine the exponent un by a data collapse. 

The numerical estimates of the critical exponents along the two horizontal sections 
through the mixed phase at k — 1 and k — 1/2 are shown in Fig. |3J The measured 
values are joined by solid lines which serve as a guide to the eye. The corresponding 
theoretical predictions are indicated by dashed lines. In the mixed phase LI the 
exponents S = (dn + a — no) /2a, 9 = k — 1, and z = a j ' k arise from the scaling relations 
and are expected to be exact. However, for the exponent v\\ we can only apply the 
one-loop expansion derived in the previous section. For the horizontal section at k — 1 
the predicted initial slope close to the mean field is again indicated by a dashed line. 
Since the cross section at k — 1/2 runs quite close to the mean field phase, we applied 
the one-loop expansion (JHT|) for the entire scope. Of course, with gradual departure 
from the phase boundary, the higher order corrections systematically increase the error 
of this approximation. 

Yet, most of the numerical estimates are in fair agreement with the theoretical 
predictions. Systematic deviations beyond the error bars are observed in the crossover 
regime from the mixed to the incubation-dominated phase (LI— Similar crossover 
effects were observed in previous numerical studies (see [21]) and may be traced back 
to temporal finite-size effects. 

Fig. El shows the corresponding results for the two vertical sections at a = 1 
and a = 3/2, respectively. In both cases the estimates for 5, 6 and z are again in 
fair agreement with the theoretical predictions, showing systematic deviations in the 
crossover region between the phases LI and L. As for the exponent V\\, the one-loop 
results (J51|) again yield initial slopes close to the mean field which now coincide with 
the numerical estimates within error bars. 

To summarize, the combination of spatial and temporal long-range interactions 
in a spreading process is a numerically challenging task. Despite of the considerable 
numerical effort the precision of our estimates is limited. Nevertheless, we consider the 
numerical evidence as sufficient to validate the two scaling relations inside the mixed 
phase. Moreover, the estimates seem to be consistent with the one-loop calculation 
of the previous section. Because of crossover effects in finite-time simulations it is 
impossible to verify the exact location of the phase boundaries. Nevertheless, all results 
are consistent with the generic phase diagram presented in Fig. 121 




(59) 
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Figure 5. Estimates of the critical exponents along the vertical sections at 
a = 1 (upper panel) and a = 1.5 (lower panel). 

4-5. Subcritical behaviour 

As first noticed in [23] the effect caused by long-range interactions can also be perceived 
within the absorbing phase (i.e., where p < p c and thus the infection becomes extinct). 
Instead of an exponential decay, which would be typical for normal DP, the extinction 
of the epidemic is now asymptotically characterized by a power law. In fact, in the 
absorbing phase the density of infected sites decreases so rapidly that interactions by 
reinfection become irrelevant, i.e. the process is asymptotically dominated by sites that 
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are infected once but then cease to newly infect other sites. Recalling that a single 
infection process complies with the distribution 



P(Ax, At) oc \Ax\- a - a Ar l - K |Ax|, At > 1 (60) 

it is evident that the last active sites generated before the infection dies out give rise 
to a temporal distribution which is simply given by the temporal component of (JoTlj) . 
Integrating out the spatial component of (J60|) . one obtains power laws for P a (t) and 
N(t) instead of exponential decays: 



Pa(t) ~ t 



-1-K 



Note that the corresponding subcritical exponents become independent of the control 
parameters a and p. 

Choosing p close to p c one observes a prolonged crossover from an almost critical 
behaviour to the algebraic decay described by (pTTjl . As we discussed in the previous 
section, for N(t) this crossover is characterized by a universal scaling function that 
involves un. We validated these findings numerically at several points within the mixed 
phase (LI in Fig. EJ). 

Interestingly, this type of asymptotic algebraic decay in the subcritical regime can 
be observed even within the DP phase (cf. Fig. |2j). This demonstrates again that the 
operators V CT and <9f differ from their short-range counterparts V 2 and dt for any finite 
a and k as they retain their long-range character. We note that this algebraic decay in 
the subcritical phase is not in contradiction with the DP hypothesis. 



5. Conclusions 



In this work we have studied the combined effect of spatial and temporal Levy flights in 
a model of non-linear interaction. Field-theoretical considerations show us that only 
the free propagator is affected by the long-range terms, while the structure of the 
Feynman-diagrams remains unaltered. Furthermore the nonlocal operators which are 
responsible for Levy flights and incubation times do not renormalize themselves, but 
instead renormalize their short-range counterparts. The calculation leads to two exact 
scaling relations which determine the three exponents 5 = (dn + a — ko~)/2ct, z = a / ' k 
and 9 — K — 1. The only non-trivial exponent measured in the simulations is u\\, which 
we calculated in the vicinity of the upper critical dimension to one-loop order. 

For a — 2, k — 1 the long-range operators <9f and V CT applied to a plane wave 
act in the same way as their short-range counterparts dt and V 2 . In a non- interacting 
model like anomalous diffusion [24], they would indeed be equivalent. However, it turns 
out that in an interacting theory, where fluctuations play an important role, long-range 
interactions are still relevant beyond this threshold. This issue is also evident from the 
location of the phase boundaries in the generic phase diagram. 
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To verify our theoretical results we have estimated the exponents along four cross 
sections through the phase diagram by numerical simulations in one spatial dimension. 
Compared with the numerical results published almost ten years ago, the present 
estimates have about the same (sometimes even less) accuracy despite of a much higher 
numerical effort. However, in comparing these results one has to keep in mind that the 
combination of spatial and temporal Levy flights is numerically much more challenging 
than the previously studied cases. Nevertheless our numerical estimates are in fair 
agreement with the predictions of the field theory. 

What can we learn from these results concerning epidemic spreading in nature? 
Recently Geisel and collaborators were able to show that individual dollar notes traced 
in space and time approximately move according to spatio-temporal Levy flights with 
exponents a ~ k ~ 1.5 [27], hence we have good reasons to expect that similar long- 
range effects will be relevant in the context of epidemic spreading. 

As a qualitative conclusion of the present study, we observe that long-range 
interactions generally decrease the critical threshold p c from where on the disease 
spreads. On the other hand, we find that the exponent f3 increases with increasing 
interaction range, meaning that the divergency of the susceptibility with respect to 
the infection rate becomes less pronounced in the presence of Levy flights, reaching a 
linear response (f3 — 1) in the mean field limit. Thus we arrive at the conclusions that 
long-range transport mechanisms of modern societies such as aviation may facilitate the 
spreading of epidemics around the globe but at the same time the onset of spreading as 
a function of the infection rate becomes less dramatic, allowing for a better containment 
of the disease. However, realistic epidemic spreading is a highly inhomogenous process, 
smearing out the phase transition over a certain range, and therefore it is still an open 
question how robust these qualitative predictions are. 
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